Stabilizacja metody elementów skończonych
Podczas wykonywania obliczeń metodą elementów skończonych wykonujemy następujące czynności
- Ustalamy obszar na którym chcemy wykonać nasze obliczenia i generujemy siatkę elementów pokrywających dany obszar, na przykład obszar \( \Omega=[0,1]^2 \) pokryty siatką 4 na 4 elementy (patchem 16 elementów na którym rozepniemy później funkcje bazowe B-spline)
- Wybieramy równanie różniczkowe w formie silnej które modeluje zjawisko fizyczne które chcemy zasymulować na naszym obszarze, na przykład \( -\Delta u=f \)
- Ustalamy warunki brzegowe które opisują co dzieje się z brzegiem naszego obszaru, na przykład \( u=0 \textrm{ dla } x=0,y\in(0,1); u=1 \textrm{ dla } x=1,y\in(0,1), \frac{\partial u}{\partial y}=0 \textrm{ dla }x\in(0,1),y\in \{ 0,1 \} \)
- Przekształcamy sformułowanie silne na tak zwane sformułowanie słabe, na przykład \( b(u,v)=l(v); b(u,v)=\int_{\Omega} \nabla u \cdot \nabla v dxdy; l(v) = \int_{\Omega} x v(x,y) dxdy \) gdzie x oznacza rozszerzenie warunku brzegowego na cały obszar (czyli taka funkcja która przyjmuje 0 na lewym brzegu obszaru i 1 na prawym, a jej pochodne w kierunku y są równe 0, tak jak wymagają tego warunki brzegowe).
- Ustalamy skończenie wymiarową bazę funkcji służących do aproksymacji naszego rozwiązania, na przykład używając wektorów węzłów na naszej siatce elementów [ 0 0 0 1 2 3 4 4 4 ] na osi x oraz [ 0 0 0 1 2 3 4 4 4 ] na osi y, co odpowiada bazie funkcji B-spline \( B^x_{i,2}B^y_{j,2}, i,j=1,...,6 \). Tworzą one przestrzeń aproksymacyjną \( U_h = span \{ B^x_{i,2}B^y_{j,2} \}_{i,j=1,...,6} \).
- Ustalamy skończenie wymiarową bazę funkcji służących do testowania sformułowania silnego, na przykład taką samą jak do aproksymacji, używając wektorów węzłów na naszej siatce elementów [ 0 0 0 1 2 3 4 4 4 ] na osi x oraz [ 0 0 0 1 2 3 4 4 4 ] na osi y, co odpowiada bazie funkcji B-spline \( B^x_{k,2}B^y_{l,2}k,l=1,...,6 \). Tworzą one przestrzeń testującą \( V_h = span \{ B^x_{k,2}B^y_{l,2} \}_{k,l=1,...,6} \).
- W ten sposób tworzymy dyskretną wersję naszego problemu \( b_h(u_h,v_h)=l(v_h) \forall v_h \in V_h \), w której zastępujemy idealne rozwiązanie \( u \) rozwiązaniem przybliżonym \( u_h = \sum_{i,j=1,...,6} u_{i,j} B^x_{i,2}(x)B^y_{j,2}(x) \), oraz jako funkcje testujące wybieramy elementy skończenie wymiarowej przestrzeni testującej \( v_h = \sum_{k,l=1,...,6} v_{k,l} B^x_{k,2}(x)B^y_{l,2}(x) \).
- Generujemy układ równań liniowych wstawiając do sformułowania słabego nasze funkcje bazowe służące do aproksymacji i testowania \( \{ \int_{\Omega} \frac{\partial B^x_{i,2} } {\partial x } B^y_{j,2} \frac{\partial B^x_{k,2} }{\partial x} B^y_{l,2}+ B^x_{i,2} \frac{\partial B^y_{j,2} }{\partial y } B^x_{k,2} \frac{\partial B^y_{l,2} }{\partial y } \}_{i,j=1,...,6;k,l=1,...,6 } \{ u_{k,l} \}_{k,l=1,...,6} = \{ \int_{\Omega} x B^x_{k,2}B^y_{l,2} \}_{k,l=1,...,6} \)
- Rozwiązujemy uzyskany układ równań liniowych, używając na przykład algorytmu LU faktoryzacji z pivotingiem
- Modyfikujemy nasze rozwiązanie odejmując od niego przesunięty warunek brzegowy \( u_h = u_h-x \).
W rozdziale tym skoncentrujemy się na zagadnieniu stabilności metody elementów skończonych. Najważniejsze twierdzenie matematyczne pozwalające badać stabilność metody elementów skończonych zostało niezależnie zaproponowane przez prof. Ivo Babuśke, prof. Franco Brezzi, oraz prof. Olga Ładyżenska [1] [2].
Odkryli oni mniej więcej w tym samym czasie, niezależnie od siebie, równoważne warunki, zwane obecnie warunkiem inf-sup ("inf-sup condition"). Warunek inf-sup jest do dzisiaj używany przez naukowców do badania stabilności metody elementów skończonych. Warunek ten zdefiniować można w abstrakcyjnych nieskończenie wymiarowych przestrzeniach matematycznych
\( \inf_{u \in U } \sup_{v \in V } \frac{b_h(u,v)}{\|u\| \|v\| } = \gamma > 0 \)
lub w przestrzeniach skończenie wymiarowych, wynikających z naszego przybliżenia rozwiązania za pomocą bazy funkcji B-spline, używanej do aproksymacji i testowani rozwiązania
\( \inf_{u \in U } \sup_{v \in V } \frac{b_h(u,v)}{\|u\| \|v\| } = \gamma > 0 \)
Może zdażyć się że nawet jeśli matematycy udowodnili spełnienie warunku inf-sup na przestrzeniach nieskończenie wymiarowych, to w momencie gdy zaczniemy rozwiązywać nasz problem w przestrzeniach skończenie wymiarowych, to warunek inf-sup nie będzie już spełniony. Na przykład, jeśli wybierzemy taką skończenie wymiarową przestrzeń testującą, że supremum z warunku inf-sup nie będzie spełnione,
\( \inf_{u \in U } \sup_{v_h \in V_h } \frac{b_h(u,v_h)}{\|u\| \|v_h\| } = \gamma > 0 \)
to wówczas nawet gdyby przestrzeń aproksymacyjna była nieskończenie wymiarowa, to i tak nasza symulacja będzie niestabilna.
Problem którego dotyczy warunek inf-sup polega na tym, iż wybierając nasze skończenie wymiarowe zbiory (tak zwane bazy) funkcji do aproksymacji rozwiązania oraz do testowania naszych równań, nie mamy pewności czy rozwiązanie uzyskane za ich pomocą na komputerze będzie dostatecznie dokładne. Żeby to sprawdzić, należy zastosować bardzo zaawansowane narzędzia matematyczne, takie jak warunek inf-sup. Dla dużej klasy problemów rozwiązywanych na komputerach za pomocą metody elementów skończonych warunek ten jest spełniony, i uzyskujemy całkiem dokładne rozwiązania. Ale jest też dość duża klasa problemów dla których warunek inf-sup nie jest spełniony.
Warunek ten napisać można w sposób następujący. Wybieramy wersję sformułowaną przez prof. Ivo Babuśkę (wszystkie wersje tego warunku opisane są w pracy Leszka Demkowicza, ICES Report 0608, 2006 "Babuśka <=> Brezzi?" ) [3]:
\( \inf_{u_h \in U_h } \sup_{v_h \in V_h } \frac{b_h(u_h,v_h)}{\|u_h\| \|v_h\| } = \gamma_h > 0 \)
Jeśli warunek ten jest spełniony, będzie dało się rozwiązać układ równań wygenerowany za pomocą metody elementów skończonych i rozwiązanie to będzie całkiem dobrym przybliżeniem idealnego rozwiązania \( u \approx u_h = \sum_{i,j=1,...,6} u_{i,j} B^x_{i,2}(x)B^y_{j,2 }(x) \).
Jeśli warunek ten nie jest spełniony, wówczas albo wygenerowany układ równań nie będzie możliwy do rozwiązania (na przykład algorytm eliminacji Gaussa przewróci się ponieważ na przekątnej macierzy powstanie numeryczne zero (bardzo mała liczba przez którą nie będzie się dało podzielić wiersza w macierzy), albo uzyskane rozwiązanie układu równań będzie niepoprawnym rozwiązaniem naszego problemu, np. powstaną oscylacje numeryczne).
W takiej sytuacji konieczna jest stabilizacja naszego problemu, która polega na
- modyfikacji funkcjonału \( b_h(u_h,v_h) \) (na dodaniu do niego jakichś dodatkowych członów stabilizujących, które w przypadku przestrzni aproksymacyjnej i testującej nieskończenie wymiarowych, w idealnym przypadku, są równe zero, ale w przypadku przestrzni skończenie wymiarowych nie są równe zero i poprawiają stabilność naszego problemu - przykładem tej metody jest stabilizacja SUPG dla równań adwekcji-dyfuzji), lub
- na modyfikacji równań problemu i normy w której próbujemy spełnić warunek inf-sup \( U_h \) (na przykład w metodzie Discontinuous Galerkin DG łamiemy funkcje i modyfikujemy równanie dodając człony równe zero w przestrzeni nieskończenie wymiarowej ale nie równe zero w przestrzeni aproksymacyjnej, modyfikujemy równanie na poziomie dyskretnym i normę co daje spełnienie warunku inf-sup również na poziomie przestrzeni dyskretnych skończenie wymiarowych), lub
- na takim przeformułowaniu problemu żeby przestrzeń aproksymacyjna i testująca były inne i żeby dało się zwiększać rozmiar przestrzeni testującej \( V_h \) (na wzięciu innej przestrzeni testującej niż przestrzeni aproksymacyjnej, dzieki czemu warunek inf-sup będzie lepiej przybliżony - tak działa na przykład metoda minimalizacji reziduum).
Wszystkie te metody stabilizacji umożliwiają matematyczne udowodnienie spełnienia warunku inf-sup. Dowody spełniania tego warunku są dość zaawansowane matematycznie, tak więc ograniczymy się jedynie do podania kilku metod stabilizacji (SUPG, minimalizacji reziduum oraz metoda DG) które gwarantują spełnianie tego warunku, ale nie podamy wyjaśnienia dlaczego ten warunek jest spełniony, ponieważ wymagałoby to wprowadzenia zaawansowanych narzędzi matematycznych.
Metoda minimalizacji reziduum leży również u podstaw metody Discontinuous Petrov-Galerkin (DPG) cieszącej się dużym powodzeniem metody stabilizacji adaptacyjnych obliczeń elementami skończonymi, zaproponowanej przez prof. Leszka Demkowicza. W metodzie DPG przestrzenie aproksymacyjne i testujące są połamane, w celu umożliwienia lokalnej statycznej kondensacji na poszczególnych elementach skończonych [4] [5] [6].
Bibliografia
1. Ivo Babuśka: Error-Bounds for Finite Element Method, Numerische Mathematik 1971, dostęp:18.10.20192. Franco Brezzi: On the Existence, Uniqueness and Approximation of Saddle-Point Problems Arising form Lagrange Multipliers, R.A.I.R.O., Italy 1974, dostęp:18.10.2019
3. Leszek Demkowicz: Babuśka <=> Brezzi, The University of Texas at Austin, USA 2006, dostęp:18.10.2019
4. T.J.R. Hughes, L.P. Franca, M. Mallet: A new finite element formulation for fluid dynamics: VI. Convergence analysis of the generalized SUPG formulation for linear time-dependent multidimensional advective–diffusive systems, Computer Methods in Applied Mechanics and Engineering, Elsevier 1987, dostęp:18.10.2019
5. Daniele Antoni Di Pietro, Alexandre Ern: Mathematical Aspects of Discontinuous Galerkin Methods, Springer, Paris 2011, dostęp:18.10.2019
6. Leszek Demkowicz, Jay Gopalakrishnan: Discontinuous Petrov–Galerkin (DPG) Method. Part 2. Fundamentals, Encyclopedia of Computational Mechanics, Wiley 2017, dostęp:18.10.2019